package EvaTCPL;

public class EvaTCPLndist extends EvaTCPLaexpression
{

  public EvaTCPLndist(EvaTCPLaexpression e1, EvaTCPLaexpression e2)
    {
      super(e1,e2);
    }

  /** Evaluate this arithmetic expression. */
  public double evaluate()
    {
      return normalInverse(Math.random())*sons[1].evaluate()+sons[0].evaluate(); 
    }

  //      ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3, 477-
  //      484.
  //
  //      Produces the normal deviate Z corresponding to a given lower
  //      tail area of P; Z is accurate to about 1 part in 10**7.
  //
  //      The hash sums below are the sums of the mantissas of the
  //      coefficients.   They are included for use in checking
  //      transcription.
  
  private static final double ZERO = 0.0, ONE = 1.0, HALF = 0.5,
      SPLIT1 = 0.425, SPLIT2 = 5.0,
      CONST1 = 0.180625, CONST2 = 1.6;

  //      Coefficients for P close to 0.5
  private static final double A0 = 3.3871327179E+00, A1 = 5.0434271938E+01,
    A2 = 1.5929113202E+02, A3 = 5.9109374720E+01,
    B1 = 1.7895169469E+01, B2 = 7.8757757664E+01,
    B3 = 6.7187563600E+01;
  //                          HASH SUM AB    32.3184577772

  //      Coefficients for P not close to 0, 0.5 or 1.
  private static final double C0 = 1.4234372777E+00, C1 = 2.7568153900E+00,
    C2 = 1.3067284816E+00, C3 = 1.7023821103E-01,
    D1 = 7.3700164250E-01, D2 = 1.2021132975E-01;
  //                          HASH SUM CD  15.7614929821

  //      Coefficients for P near 0 or 1.
  private static final double E0 = 6.6579051150E+00, E1 = 3.0812263860E+00,
    E2 = 4.2868294337E-01, E3 = 1.7337203997E-02,
    F1 = 2.4197894225E-01, F2 = 1.2258202635E-02;
  //                          HASH SUM EF    19.4052910204

  /** Inverse of N(0,1) distributed normal distribution.
      @param P The probability (between 0 and 1).
  */
  private static double normalInverse(double P)// throws ArithmeticException
    {
      double   PPND7, Q, R;

      Q = P - HALF;
      if (Math.abs(Q)<=SPLIT1) {
	R = CONST1 - Q * Q;
	return Q * (((A3 * R + A2) * R + A1) * R + A0) /
	  (((B3 * R + B2) * R + B1) * R + ONE);
      }
      else if (Q<ZERO)
        R = P;
      else
        R = ONE - P;

      if ( R == ZERO )
        R = Double.MIN_VALUE;       // Constrain R to greater than zero.
      else if (R<ZERO)
	System.out.println("P="+P);
      //        throw new ArithmeticException("P="+P);

      R = Math.sqrt( -Math.log(R) );

      if (R <= SPLIT2){
        R = R - CONST2;
        PPND7 = (((C3 * R + C2) * R + C1) * R + C0) /
	  ((D2 * R + D1) * R + ONE);
      }
      else {
        R = R - SPLIT2;
        PPND7 = (((E3 * R + E2) * R + E1) * R + E0) /
	  ((F2 * R + F1) * R + ONE);
      }

      if (Q<ZERO)
        PPND7 = - PPND7;

      return PPND7;
    }

  public String toString()
    {
      return "ndist("+sons[0].toString()+", "+sons[1].toString()+")";
    }


}
